Logistic and Autologistic in R
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(RColorBrewer)
sl <- read_sf("./SriLankaNYCFinal/SriLankaCorrectedFinal_UTM.shp")
#create a poor Secretariat on not variable
sl$POOR[is.na(sl$POOR)] <- 0
sl$POORFac <- factor(sl$POOR)
#ESDA
#map poor or not
plot(sl["POOR"],main='Poverty by Secretariat (binary)', breaks="jenks")
#map percentage poor
plot(sl["PCTPOOHH"],main='Poverty by Secretariat', breaks="jenks")
#histograms to see the distribution
hist(sl$POOR)
hist(sl$PCTPOOHH, breaks=20)
summary(sl)
## OBJECTID ID_0 ISO NAME_0
## Min. :36745 Min. :217 Length:236 Length:236
## 1st Qu.:36822 1st Qu.:217 Class :character Class :character
## Median :36898 Median :217 Mode :character Mode :character
## Mean :36895 Mean :217
## 3rd Qu.:36967 3rd Qu.:217
## Max. :37033 Max. :217
## ID_1 NAME_1 ID_2 NAME_2
## Min. : 2.00 Length:236 Min. : 20.00 Length:236
## 1st Qu.: 6.00 Class :character 1st Qu.: 96.75 Class :character
## Median :12.00 Mode :character Median :172.50 Mode :character
## Mean :12.13 Mean :169.53
## 3rd Qu.:17.00 3rd Qu.:242.25
## Max. :23.00 Max. :308.00
## HASC_2 CCN_2 CCA_2 TYPE_2
## Length:236 Min. :0 Length:236 Length:236
## Class :character 1st Qu.:0 Class :character Class :character
## Mode :character Median :0 Mode :character Mode :character
## Mean :0
## 3rd Qu.:0
## Max. :0
## ENGTYPE_2 NL_NAME_2 VARNAME_2 Shape_Leng
## Length:236 Length:236 Length:236 Min. :0.2102
## Class :character Class :character Class :character 1st Qu.:0.4752
## Mode :character Mode :character Mode :character Median :0.6094
## Mean :0.6897
## 3rd Qu.:0.7957
## Max. :2.4245
## Shape_Area FID_1 ADM3CODE ADM1CODE
## Min. :0.001370 Min. : 0.00 Min. : 1.00 Min. :1.000
## 1st Qu.:0.006363 1st Qu.: 60.75 1st Qu.: 61.75 1st Qu.:2.000
## Median :0.010850 Median :196.50 Median :197.50 Median :6.000
## Mean :0.015412 Mean :161.38 Mean :162.38 Mean :4.737
## 3rd Qu.:0.018074 3rd Qu.:261.25 3rd Qu.:262.25 3rd Qu.:7.000
## Max. :0.087908 Max. :322.00 Max. :323.00 Max. :9.000
## ADM1NAME ADM2NAME ADM2CODE ADM3NAME
## Length:236 Length:236 Min. :1103 Length:236
## Class :character Class :character 1st Qu.:2208 Class :character
## Mode :character Mode :character Median :6110 Mode :character
## Mean :4922
## 3rd Qu.:7175
## Max. :9233
## NUMHH NUMPER NUMPER_U NUMPER_R
## Min. : 2990 Min. : 12450 Min. : 0 Min. : 0
## 1st Qu.: 8242 1st Qu.: 32880 1st Qu.: 0 1st Qu.: 29711
## Median :12660 Median : 51391 Median : 0 Median : 45347
## Mean :15745 Mean : 65833 Mean : 9299 Mean : 52848
## 3rd Qu.:18097 3rd Qu.: 74626 3rd Qu.: 0 3rd Qu.: 61871
## Max. :86581 Max. :377396 Max. :377396 Max. :209741
## NUMPER_EST NUMPER_M NUMPER_F NUMPER17
## Min. : 0.0 Min. : 6508 Min. : 5942 Min. : 4184
## 1st Qu.: 0.0 1st Qu.: 16170 1st Qu.: 16622 1st Qu.: 11800
## Median : 92.5 Median : 25566 Median : 26697 Median : 17878
## Mean : 3685.6 Mean : 32542 Mean : 33291 Mean : 21475
## 3rd Qu.: 2173.0 3rd Qu.: 36353 3rd Qu.: 37754 3rd Qu.: 24637
## Max. :142483.0 Max. :203606 Max. :173790 Max. :113606
## ETSINPOP ETSLTPOP ETHINTPO ETHSLMPO
## Min. : 9026 Min. : 0.0 Min. : 0 Min. : 0
## 1st Qu.: 28267 1st Qu.: 137.8 1st Qu.: 15 1st Qu.: 76
## Median : 44298 Median : 681.0 Median : 284 Median : 1344
## Mean : 54961 Mean : 2477.9 Mean : 3484 Mean : 4464
## 3rd Qu.: 63581 3rd Qu.: 1679.8 3rd Qu.: 2070 3rd Qu.: 4850
## Max. :206081 Max. :117510.0 Max. :141066 Max. :120109
## ETHBURPO ETHMALPO ETHOTHPO RELBUDPO
## Min. : 0.00 Min. : 0.00 Min. : 0.0 Min. : 4463
## 1st Qu.: 3.00 1st Qu.: 1.75 1st Qu.: 2.0 1st Qu.: 27427
## Median : 11.50 Median : 9.00 Median : 11.5 Median : 42501
## Mean : 137.16 Mean : 192.06 Mean : 117.4 Mean : 51381
## 3rd Qu.: 49.75 3rd Qu.: 47.00 3rd Qu.: 50.5 3rd Qu.: 61387
## Max. :3380.00 Max. :8680.00 Max. :5667.0 Max. :195882
## RELHINPO RELISLPO RELCATPO RELOTHPO
## Min. : 0 Min. : 1.0 Min. : 11.0 Min. : 0.00
## 1st Qu.: 119 1st Qu.: 115.5 1st Qu.: 211.5 1st Qu.: 1.00
## Median : 950 Median : 1432.0 Median : 620.5 Median : 7.50
## Mean : 4966 Mean : 4764.7 Mean : 4685.6 Mean : 35.06
## 3rd Qu.: 3735 3rd Qu.: 5223.0 3rd Qu.: 2064.5 3rd Qu.: 25.00
## Max. :140176 Max. :134271.0 Max. :98254.0 Max. :1217.00
## NOHOSU20 NOCLQ200 NOINS200 NONHU200
## Min. : 3651 Min. : 0.0 Min. : 14.00 Min. : 352
## 1st Qu.: 9311 1st Qu.: 13.0 1st Qu.: 47.00 1st Qu.: 1123
## Median :14185 Median : 28.5 Median : 70.50 Median : 1763
## Mean :17071 Mean : 144.0 Mean : 85.57 Mean : 2276
## 3rd Qu.:19943 3rd Qu.: 67.5 3rd Qu.:109.25 3rd Qu.: 2610
## Max. :66189 Max. :9029.0 Max. :628.00 Max. :27385
## TOTBN200 AGHOLB40 AGHLDA40 TOTAGHOL
## Min. : 4173 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.:10651 1st Qu.: 1087 1st Qu.: 4240 1st Qu.: 7054
## Median :15949 Median : 2976 Median : 6004 Median :10300
## Mean :19577 Mean : 4940 Mean : 6500 Mean :11440
## 3rd Qu.:22854 3rd Qu.: 6127 3rd Qu.: 8441 3rd Qu.:13769
## Max. :95311 Max. :32141 Max. :22477 Max. :42982
## NOPNOAL2 NOPOOHG2 NOPOHGOL NOPOOLO2
## Min. : 0.0 Min. : 0 Min. : 0 Min. : 0.0
## 1st Qu.: 152.2 1st Qu.: 687 1st Qu.: 1002 1st Qu.: 472.5
## Median : 312.0 Median :1248 Median : 1888 Median : 1286.5
## Mean : 468.1 Mean :1505 Mean : 2247 Mean : 2126.9
## 3rd Qu.: 612.2 3rd Qu.:2146 3rd Qu.: 3169 3rd Qu.: 3174.8
## Max. :5443.0 Max. :8418 Max. :10052 Max. :11436.0
## NOPTOT20 NAHLT1AC AAHLT1AC NAHB1A2A ACAHB1A2
## Min. : 0 Min. : 0 Min. : 0.0 Min. : 0 Min. : 0
## 1st Qu.: 4173 1st Qu.:1111 1st Qu.: 598.0 1st Qu.:1184 1st Qu.:1524
## Median : 6016 Median :1884 Median : 982.5 Median :1810 Median :2393
## Mean : 6347 Mean :2144 Mean :1116.7 Mean :2007 Mean :2605
## 3rd Qu.: 8168 3rd Qu.:2926 3rd Qu.:1578.5 3rd Qu.:2738 3rd Qu.:3536
## Max. :20992 Max. :6787 Max. :3626.0 Max. :7033 Max. :9077
## NAHA2AC2 ACAHA2AC TNAH2001 TACAH200
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 1051 1st Qu.: 4230 1st Qu.: 4220 1st Qu.: 7242
## Median : 1808 Median : 7221 Median : 6156 Median :11898
## Mean : 2328 Mean : 8734 Mean : 6472 Mean :12456
## 3rd Qu.: 3222 3rd Qu.:11911 3rd Qu.: 8408 3rd Qu.:17118
## Max. :10840 Max. :36941 Max. :21065 Max. :46643
## TEAAC200 RUBAC200 COCACIP2 COCACNIP
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 191.2 1st Qu.: 40.75
## Median : 0.0 Median : 7.0 Median : 510.0 Median : 153.00
## Mean : 929.9 Mean : 488.5 Mean : 1628.2 Mean : 306.21
## 3rd Qu.: 694.8 3rd Qu.: 245.5 3rd Qu.: 1476.2 3rd Qu.: 437.75
## Max. :11253.0 Max. :7262.0 Max. :17939.0 Max. :2599.00
## CASHAC20 TNONCA20 TNOCIM20 DCMPROD2
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 2.0 1st Qu.: 514.2 1st Qu.: 114.5 1st Qu.: 216.0
## Median : 20.5 Median : 1276.5 Median : 237.0 Median : 571.5
## Mean : 155.9 Mean : 2714.2 Mean : 446.3 Mean : 1003.2
## 3rd Qu.: 177.8 3rd Qu.: 4022.8 3rd Qu.: 587.0 3rd Qu.: 1469.0
## Max. :4369.0 Max. :20028.0 Max. :2902.0 Max. :11254.0
## TNOB2001 TNOBIM20 DBMPROD2 TNOGS200
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 163.8 1st Qu.: 4.0 1st Qu.: 7.0 1st Qu.: 87.75
## Median : 385.0 Median : 40.5 Median : 88.0 Median : 229.00
## Mean : 923.3 Mean : 154.5 Mean : 253.2 Mean : 526.47
## 3rd Qu.:1180.2 3rd Qu.: 146.0 3rd Qu.: 307.8 3rd Qu.: 663.00
## Max. :6389.0 Max. :1675.0 Max. :2730.0 Max. :3654.00
## TNOS2001 TNOC2001 SAM1000R SAM700R2
## Min. : 0.00 Min. : 0 Min. : 0.00 Min. : 180
## 1st Qu.: 3.75 1st Qu.: 3347 1st Qu.: 4.00 1st Qu.: 1490
## Median : 44.50 Median : 7886 Median : 17.00 Median : 2864
## Mean : 229.68 Mean : 32453 Mean : 34.59 Mean : 3320
## 3rd Qu.: 234.25 3rd Qu.: 26224 3rd Qu.: 44.25 3rd Qu.: 4588
## Max. :3631.00 Max. :606821 Max. :379.00 Max. :11675
## SAM500R2 SAM400R2 SAM350R2 SAM250R2
## Min. : 0.00 Min. : 0 Min. : 202 Min. : 111.0
## 1st Qu.: 2.00 1st Qu.: 3 1st Qu.: 695 1st Qu.: 454.0
## Median : 11.00 Median : 13 Median :1080 Median : 772.5
## Mean : 20.29 Mean :1275 Mean :1404 Mean : 889.3
## 3rd Qu.: 29.00 3rd Qu.:2534 3rd Qu.:1682 3rd Qu.:1157.5
## Max. :220.00 Max. :9900 Max. :7077 Max. :4287.0
## SAM140R2 SAMTNOF2 NOPAPAR3 ASSWMAJ3
## Min. : 0.00 Min. : 1635 Min. : 0 Min. : 0
## 1st Qu.: 0.00 1st Qu.: 4578 1st Qu.: 2528 1st Qu.: 0
## Median : 0.00 Median : 6304 Median : 4508 Median : 0
## Mean : 31.03 Mean : 6961 Mean : 5422 Mean : 1948
## 3rd Qu.: 14.25 3rd Qu.: 8683 3rd Qu.: 7922 3rd Qu.: 1413
## Max. :1535.00 Max. :19405 Max. :21754 Max. :31780
## ASSWMIN3 ASSWRF3 ASSWTOT3 PHAMJI00
## Min. : 0.0 Min. : 0.0 Min. : 0 Min. : 0.0
## 1st Qu.: 284.8 1st Qu.: 140.5 1st Qu.: 1719 1st Qu.: 0.0
## Median : 780.0 Median : 746.5 Median : 3256 Median : 0.0
## Mean : 1430.0 Mean :1302.7 Mean : 4689 Mean : 1401.5
## 3rd Qu.: 1763.2 3rd Qu.:2010.8 3rd Qu.: 5613 3rd Qu.: 587.2
## Max. :12545.0 Max. :6804.0 Max. :31780 Max. :29104.0
## PHAMII00 PHARF00 PHATO00 PHAMJI01
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0
## 1st Qu.: 106.5 1st Qu.: 29.75 1st Qu.: 926.5 1st Qu.: 0
## Median : 427.5 Median : 422.50 Median : 2005.0 Median : 0
## Mean : 878.4 Mean : 972.22 Mean : 3252.2 Mean : 2797
## 3rd Qu.: 1096.0 3rd Qu.:1623.50 3rd Qu.: 4090.8 3rd Qu.: 1075
## Max. :12469.0 Max. :5417.00 Max. :29104.0 Max. :59047
## PHAMII01 PHARF01 PHATO01 DISROADS
## Min. : 0 Min. : 0 Min. : 0 Min. : 0.050
## 1st Qu.: 216 1st Qu.: 37 1st Qu.: 1559 1st Qu.: 0.620
## Median : 828 Median : 577 Median : 3474 Median : 1.050
## Mean : 1419 Mean :1423 Mean : 5640 Mean : 1.521
## 3rd Qu.: 1919 3rd Qu.:2293 3rd Qu.: 6873 3rd Qu.: 1.900
## Max. :13522 Max. :8844 Max. :59047 Max. :11.940
## DISTOWN MAHARF YALARF TOTRF
## Min. : 0.00 Min. :143.0 Min. : 727.0 Min. :1178
## 1st Qu.:10.00 1st Qu.:376.0 1st Qu.: 973.8 1st Qu.:1739
## Median :18.00 Median :558.0 Median :1049.0 Median :1960
## Mean :19.42 Mean :564.7 Mean :1073.4 Mean :2044
## 3rd Qu.:27.00 3rd Qu.:745.0 3rd Qu.:1195.0 3rd Qu.:2402
## Max. :48.00 Max. :968.0 Max. :1405.0 Max. :2979
## AVE_ELEV MAX_ELEV AVE_SLOP MAX_SLOP
## Min. : 0.155 Min. : 5.878 Min. : 2.962 Min. : 0.976
## 1st Qu.: 1.171 1st Qu.: 126.250 1st Qu.: 37.417 1st Qu.:21.000
## Median : 3.397 Median : 365.500 Median : 101.088 Median :35.000
## Mean : 5.561 Mean : 616.380 Mean : 230.341 Mean :39.022
## 3rd Qu.: 9.343 3rd Qu.: 950.250 3rd Qu.: 287.656 3rd Qu.:53.500
## Max. :18.785 Max. :2503.000 Max. :1670.571 Max. :85.000
## PCTPOOHH POVGROUP NUMPOHH FGT_0
## Min. : 1.001 Min. :1.000 Min. : 367 Min. : 1.143
## 1st Qu.:23.370 1st Qu.:3.000 1st Qu.: 2401 1st Qu.:28.167
## Median :29.555 Median :4.000 Median : 3386 Median :34.892
## Mean :27.737 Mean :3.962 Mean : 3698 Mean :32.614
## 3rd Qu.:32.803 3rd Qu.:5.000 3rd Qu.: 4707 3rd Qu.:38.345
## Max. :45.650 Max. :6.000 Max. :10215 Max. :52.854
## POORPER PCTSAMHH SAMBENGP SAMMISAL
## Min. : 1531 Min. : 5.197 Min. :1.000 Min. :1.000
## 1st Qu.: 9659 1st Qu.: 37.708 1st Qu.:5.000 1st Qu.:1.000
## Median :14187 Median : 48.712 Median :6.000 Median :3.000
## Mean :15406 Mean : 48.371 Mean :5.475 Mean :2.581
## 3rd Qu.:20407 3rd Qu.: 58.604 3rd Qu.:6.000 3rd Qu.:4.000
## Max. :45218 Max. :100.525 Max. :6.000 Max. :5.000
## POPDENS id2 SMHOLDPAG AGOPWOLA
## Min. : 8.61 Min. : 1.00 Min. : 0.00 Min. :0.00000
## 1st Qu.: 68.83 1st Qu.: 61.75 1st Qu.:18.48 1st Qu.:0.03330
## Median :131.28 Median :197.50 Median :37.65 Median :0.05320
## Mean :168.03 Mean :162.38 Mean :35.50 Mean :0.07213
## 3rd Qu.:210.91 3rd Qu.:262.25 3rd Qu.:48.97 3rd Qu.:0.09084
## Max. :778.22 Max. :323.00 Max. :88.55 Max. :0.57604
## PACSMALL NUMAGHH POOR geometry POORFac
## Min. : 0.000 Min. :0.0000 Min. :0.0000 MULTIPOLYGON :236 0:179
## 1st Qu.: 4.524 1st Qu.:0.3197 1st Qu.:0.0000 epsg:4326 : 0 1: 57
## Median :11.544 Median :0.5433 Median :0.0000 +proj=long...: 0
## Mean :11.559 Mean :0.5379 Mean :0.2415
## 3rd Qu.:16.183 3rd Qu.:0.7564 3rd Qu.:0.0000
## Max. :42.813 Max. :1.2079 Max. :1.0000
library(spatialEco)
## Warning: package 'spatialEco' was built under R version 4.4.2
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
#run a logistic regression
sl_lmodel <- logistic.regression(sl, y="POOR",
x=c("NUMAGHH","DISTOWN","DISROADS", "ASSWMAJ3"))
sl_lmodel$model
## Logistic Regression Model
##
## rms::lrm(formula = form, data = ldata, x = TRUE, y = TRUE, penalty = bf$penalty)
##
##
## Penalty factors
##
## simple nonlinear interaction nonlinear.interaction
## 1 1 1 1
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 236 LR chi2 32.74 R2 0.189 C 0.744
## 0 179 d.f. 3.851 R2(4,236)0.112 Dxy 0.488
## 1 57 Pr(> chi2)<0.0001 R2(4,129.7)0.194 gamma 0.488
## max |deriv| 6e-11 Penalty 0.81 Brier 0.160 tau-a 0.179
##
## Coef S.E. Wald Z Pr(>|Z|) Penalty Scale
## Intercept -3.1991 0.4818 -6.64 <0.0001 0.0000
## NUMAGHH 3.0747 0.8022 3.83 0.0001 0.2773
## DISTOWN 0.0077 0.0176 0.44 0.6607 11.0418
## DISROADS 0.0988 0.1149 0.86 0.3899 1.5223
## ASSWMAJ3 -0.0001 0.0000 -1.27 0.2038 4331.3029
sl_lmodel$diagTable
## Names Value
## 1 Obs 2.360000e+02
## 2 Max Deriv 6.355094e-11
## 3 Model L.R. 3.273924e+01
## 4 d.f. 3.851091e+00
## 5 P 1.124815e-06
## 6 C 7.438008e-01
## 7 Dxy 4.876017e-01
## 8 Gamma 4.876017e-01
## 9 Tau-a 1.794086e-01
## 10 R2 1.891640e-01
## 11 R2(236) 1.265524e-01
## 12 R2(4,236) 1.116221e-01
## 13 R2(129.7) 2.182375e-01
## 14 R2(4,129.7) 1.937519e-01
## 15 Brier 1.598022e-01
## 16 g 1.061064e+00
## 17 gr 2.889443e+00
## 18 gp 1.731764e-01
## 19 PEN 1.000000e+00
## 20 AIC 2.503706e+01
sl_lmodel$coefTable
## Variable Coef StdError Wald Prob
## 1 Intercept -3.199070e+00 4.817686e-01 -6.6402618 3.131265e-11
## 2 NUMAGHH 3.074701e+00 8.022077e-01 3.8327994 1.266933e-04
## 3 DISTOWN 7.747269e-03 1.764783e-02 0.4389926 6.606669e-01
## 4 DISROADS 9.879432e-02 1.149143e-01 0.8597215 3.899426e-01
## 5 ASSWMAJ3 -5.152669e-05 4.054994e-05 -1.2706970 2.038365e-01
#save residuals and predicted probabilities
VarNamesLogModel <- rownames(sl_lmodel$model$var)[-1]
sl$LogResiduals <- sl_lmodel$Residuals[,1]
sl$LogStdResiduals <- sl_lmodel$Residuals[,2]
sl$LogProbs <- predict(sl_lmodel$model,
sf::st_drop_geometry(sl[,VarNamesLogModel]),
type='fitted')
#plot the probability of being poor predicted by autologistic model
plot(sl["LogProbs"],main='Poverty probability predicted by Secretariat
Logistic Model', breaks="jenks")
#make a dataframe of the variables in the layer
sl_coord <- st_coordinates(st_centroid(sl))
## Warning: st_centroid assumes attributes are constant over geometries
sl_df <- data.frame(sl)
sl_autolmodel <- logistic.regression(sl_df,
y="POOR",
x=c("NUMAGHH","DISTOWN","DISROADS", "ASSWMAJ3"),
autologistic=TRUE,
coords=sl_coord)
sl_autolmodel$model
## Logistic Regression Model
##
## rms::lrm(formula = form, data = ldata, x = TRUE, y = TRUE, penalty = bf$penalty)
##
##
## Penalty factors
##
## simple nonlinear interaction nonlinear.interaction
## 0.9 0.9 0.9 0.9
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 236 LR chi2 121.14 R2 0.587 C 0.929
## 0 179 d.f. 4.698 R2(5,236)0.379 Dxy 0.859
## 1 57 Pr(> chi2)<0.0001 R2(5,129.7)0.580 gamma 0.859
## max |deriv| 3e-10 Penalty 3.52 Brier 0.090 tau-a 0.316
##
## Coef S.E. Wald Z Pr(>|Z|) Penalty Scale
## Intercept -4.8686 0.7517 -6.48 <0.0001 0.0000
## NUMAGHH 3.8686 1.1322 3.42 0.0006 0.2631
## DISTOWN -0.0138 0.0243 -0.57 0.5699 10.4752
## DISROADS 0.0997 0.1622 0.61 0.5387 1.4442
## ASSWMAJ3 -0.0001 0.0001 -0.95 0.3428 4109.0347
## AutoCov 4.9305 0.6500 7.59 <0.0001 0.3142
sl_autolmodel$diagTable
## Names Value
## 1 Obs 2.360000e+02
## 2 Max Deriv 3.224159e-10
## 3 Model L.R. 1.211375e+02
## 4 d.f. 4.698289e+00
## 5 P 0.000000e+00
## 6 C 9.294325e-01
## 7 Dxy 8.588650e-01
## 8 Gamma 8.588650e-01
## 9 Tau-a 3.160115e-01
## 10 R2 5.866516e-01
## 11 R2(236) 3.924751e-01
## 12 R2(5,236) 3.794665e-01
## 13 R2(129.7) 5.961928e-01
## 14 R2(5,129.7) 5.803217e-01
## 15 Brier 8.995618e-02
## 16 g 2.386243e+00
## 17 gr 1.087257e+01
## 18 gp 2.998143e-01
## 19 PEN 9.000000e-01
## 20 AIC 1.117409e+02
sl_autolmodel$coefTable
## Variable Coef StdError Wald Prob
## 1 Intercept -4.868592e+00 7.516760e-01 -6.4769823 9.357515e-11
## 2 NUMAGHH 3.868608e+00 1.132242e+00 3.4167693 6.336895e-04
## 3 DISTOWN -1.382178e-02 2.432633e-02 -0.5681819 5.699115e-01
## 4 DISROADS 9.972160e-02 1.622130e-01 0.6147570 5.387152e-01
## 5 ASSWMAJ3 -5.237384e-05 5.520869e-05 -0.9486520 3.427976e-01
## 6 AutoCov 4.930545e+00 6.499581e-01 7.5859426 3.300784e-14
#add the residual etc to the dataframe
sl$AutoCov <- sl_autolmodel$AutoCov
sl$AutoLogResiduals <- sl_autolmodel$Residuals[,1]
sl$AutoLogStdResiduals <- sl_autolmodel$Residuals[,2]
sl$AutoLogProbs <- predict(sl_autolmodel$model,
type='fitted')
#plot the probability of being poor predicted by autologistic model
plot(sl["AutoLogProbs"],main='Poverty probability predicted by Secretariat
Autologistic Model', breaks="jenks")
OLS and GWR in R
library(spdep)
library(spgwr)
## Warning: package 'spgwr' was built under R version 4.4.3
## Loading required package: sp
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
library(RColorBrewer)
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
# OLS model
sl_pov_olsmodel <- lm(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN + ASSWMAJ3, data=sl)
summary(sl_pov_olsmodel)
##
## Call:
## lm(formula = PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN +
## ASSWMAJ3, data = sl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5146 -4.3922 -0.1625 3.1168 17.1759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.6716175 1.9800545 12.460 < 2e-16 ***
## NUMAGHH 17.5793199 1.9896718 8.835 2.62e-16 ***
## MAHARF -0.0111072 0.0024187 -4.592 7.22e-06 ***
## DISROADS -0.2471006 0.3511222 -0.704 0.4823
## DISTOWN 0.0389303 0.0455097 0.855 0.3932
## ASSWMAJ3 -0.0002558 0.0001085 -2.358 0.0192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.142 on 230 degrees of freedom
## Multiple R-squared: 0.487, Adjusted R-squared: 0.4758
## F-statistic: 43.66 on 5 and 230 DF, p-value: < 2.2e-16
AIC(sl_pov_olsmodel)
## [1] 1534.377
# GWR models
adaptive_bw <- gwr.sel(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS +
DISTOWN + ASSWMAJ3,
data=sl, method = "aic", adapt = TRUE, coords=sl_coord)
## Bandwidth: 0.381966 AIC: 1480.366
## Bandwidth: 0.618034 AIC: 1502.147
## Bandwidth: 0.236068 AIC: 1450.837
## Bandwidth: 0.145898 AIC: 1411.916
## Bandwidth: 0.09016994 AIC: 1370.652
## Bandwidth: 0.05572809 AIC: 1332.262
## Bandwidth: 0.03444185 AIC: 1321.457
## Bandwidth: 0.02178569 AIC: 1347.974
## Bandwidth: 0.04177523 AIC: 1318.577
## Bandwidth: 0.04115187 AIC: 1319
## Bandwidth: 0.04710475 AIC: 1321.653
## Bandwidth: 0.04381093 AIC: 1318.699
## Bandwidth: 0.04268552 AIC: 1318.296
## Bandwidth: 0.04272621 AIC: 1318.304
## Bandwidth: 0.04252023 AIC: 1318.269
## Bandwidth: 0.04223567 AIC: 1318.319
## Bandwidth: 0.04241153 AIC: 1318.255
## Bandwidth: 0.04234436 AIC: 1318.265
## Bandwidth: 0.04245222 AIC: 1318.26
## Bandwidth: 0.04241153 AIC: 1318.255
adaptive_bw
## [1] 0.04241153
# Optimal adaptive bandwidth for the above case is: approximately 0.04
# each regression will expand or contract to include roughly 4% of sample
fixed_bw <- ggwr.sel(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS +
DISTOWN + ASSWMAJ3,
data=sl, adapt = FALSE, coords=sl_coord)
## Bandwidth: 1.297489 CV score: 8242.629
## Bandwidth: 2.097285 CV score: 8755.026
## Bandwidth: 0.8031877 CV score: 7255.535
## Bandwidth: 0.4976927 CV score: 5854.028
## Bandwidth: 0.3088864 CV score: 4421.947
## Bandwidth: 0.1921977 CV score: 3786.395
## Bandwidth: 0.1200801 CV score: 4290.592
## Bandwidth: 0.2092028 CV score: 3803.537
## Bandwidth: 0.1950849 CV score: 3786.397
## Bandwidth: 0.1936372 CV score: 3786.24
## Bandwidth: 0.1936779 CV score: 3786.24
## Bandwidth: 0.1935965 CV score: 3786.24
## Bandwidth: 0.1936372 CV score: 3786.24
fixed_bw
## [1] 0.1936372
# 21415.7
# This means that the bandwidth for each regression using a fixed bandwidth
# will include all centroid within approx 21 km
# USING THE ADAPTIVE BISQUARE KERNEL
# Note that the input for the adapt command is the value generated
# in the calculation above for adaptive.bw
sl_adgwr_model <- gwr(PCTPOOHH ~ NUMAGHH + MAHARF + DISROADS + DISTOWN +
ASSWMAJ3, data=sl, adapt = adaptive_bw,
hatmatrix = TRUE,
coords=sl_coord)
names(sl_adgwr_model)
## [1] "SDF" "lhat" "lm" "results" "bandwidth" "adapt"
## [7] "hatmatrix" "gweight" "gTSS" "this.call" "fp.given" "timings"
summary(sl_adgwr_model)
## Length Class Mode
## SDF 236 SpatialPointsDataFrame S4
## lhat 55696 -none- numeric
## lm 11 -none- list
## results 14 -none- list
## bandwidth 236 -none- numeric
## adapt 1 -none- numeric
## hatmatrix 1 -none- logical
## gweight 1 -none- character
## gTSS 1 -none- numeric
## this.call 6 -none- call
## fp.given 1 -none- logical
## timings 12 -none- numeric
names(sl_adgwr_model$results)
## [1] "v1" "v2" "delta1" "delta2" "sigma2" "sigma2.b"
## [7] "AICb" "AICh" "AICc" "edf" "rss" "nu1"
## [13] "odelta2" "n"
sl_adgwr_model$results$AICc
## [1] 1330.992
#Note the drop in AICc compared to regression AIC
# save local coefficients and standard errors to a new dataframe
gwrmodel_results <- as.data.frame(sl_adgwr_model$SDF)
# join it with the shapefile using ID_2 to match rows for both
# join it with the spatial data
sl_gwr <- cbind(sl, as.matrix(gwrmodel_results))
sl_gwr_sf <- st_as_sf(sl_gwr)
# write your spatial data from R to a shapefile
write_sf(sl_gwr, "SriLankaGWR.shp", OVERWRITE=TRUE)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
# Repeat this section for each variable in your GWR
# first create a new variable NUMAGHH_T with significant coeff
#note that R has added a .1 after the variable to show that
#this is the coefficient and not the original variable
sl_gwr$NUMAGHH_T <- sl_gwr$NUMAGHH.1/sl_gwr$NUMAGHH_se
sl_gwr$DISROADS_T <- sl_gwr$DISROADS.1/sl_gwr$DISROADS_se
#If there are values of Max and/or Min over 1.96 we map it
#if its not significant use the OLS model coefficient
summary(sl_gwr$NUMAGHH_T)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.6777 2.3531 3.3901 4.0401 5.3722 13.5991
summary(sl_gwr$DISROADS_T)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.5373 -0.0109 0.4778 0.6146 1.4391 5.9105
#setting the coefficient so that it is only symbolized if T is significant
sl_gwr$CoeffNUMAGHH <- ifelse(abs(sl_gwr$NUMAGHH_T)>=1.96,sl_gwr$NUMAGHH_T, NA)
sl_gwr$CoeffDISROADS <- ifelse(abs(sl_gwr$DISROADS_T)>=1.96,sl_gwr$DISROADS_T, NA)
#Tmap is better for mapping multiple layers
#tmap_mode is set to view so that you can zoom in to see location names
#change the basemap from the default to Openstreetmap to see place names
#see all pallettes by running the next line without the comment #
#cols4all::c4a_palettes()
map_slr2 <- tm_shape(sl_gwr) +
tm_polygons(fill="localR2",
fill.scale = tm_scale_intervals(values="brewer.blues", style = "jenks"),
fill.legend = tm_legend(title = "R Squared", size = 0.8))+
tm_borders() + tmap_mode("view")
## ℹ tmap mode set to "view".
map_slr2